#'
#'
#' @description Based on (2.18) in the chapter2 in the Heer (2019)
#'
#'
diff1_heer <- function(y, alpha = 0.36, delta = 0.08, n = 0, beta0 = 0.96, sigma = 2){
  # local f,c1,c2,k1,k2;
  k1 <- y[1]
  c1 <- y[2]
  f <- matlab::zeros(2,1)
  f[1] <- k1^(alpha)+(1-delta)*k1-c1
  f[1] <- f[1]/(1+n)
  f[2] <- beta0^(1/sigma)*(1+alpha* ( (k1^(alpha)+(1-delta)*k1-c1) / (1+n) )^(alpha-1)-delta)^(1/sigma)*c1
  return(f)
}

